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Abstract 


The  effects  of  thermal  nonequilibrium  on  flows  about 
blunt  bodies  have  not  been  studied  numerically  in  isolation 
from  chemical  reactions.  Typically,  air  is  modeled  as  a 
perfect  gas  or  as  a  chemically  reacting  mixture.  In  the 
former  case,  significant  errors  result  at  Mach  numbers 
exceeding  about  5.  However,  below  Mach  numbers  around  8,  the 
effects  of  chemical  reactions  are  negligible.  This  study 
examines  inviscid  flow  about  a  simple  axisymmetric  blunt-body 
at  moderate  hypersonic  speeds  (Mach  numbers  between  5-8) .  A 
first-order-accurate  Roe  scheme  was  used  to  compute  solutions 
of  the  Euler  equations  assuming  different  gas  models:  perfect 
gas,  thermal  equilibrium,  and  thermal  nonequilibrium.  Two 
vibrational  relaxation  models  were  used  to  examine  thermal 
nonequilibrium.  Equilibrium  and  nonequilibrium  shock  standoff 
distances  deviated  10-15  percent  from  perfect  gas 
calculations.  Nonequilibrium  calculations  produced  striking 
temperature  and  density  gradients  in  the  stagnation  region 
with  an  isolated  "hot  spot"  near  the  shock. 
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A  COMPARISON  OF  MOLECULAR  VIBRATION  MODELING 


FOR  THERMAL  NONEQUILIBRIUM  AIRFLOW 


1.  Introduction 


High  velocities  and  Mach  numbers,  strong  shock  waves, 
high-temperature  phenomena,  and  shock/boundary-layer 
interactions  comprise  a  few  aspects  of  hypersonic  flows.  Some 
of  these  features  are  viscous  in  nature,  but  others  also  occur 
in  an  inviscid  environment.  The  inviscid  blunt-body  flow 
features  a  strong  shock  which  is  normal  to  the  nose,  a  thin 
shock  layer  and  thermo-chemical  nonequilibrium.  Viscous 
features  include  "entropy  swallowing"  by  the  boundary-layer, 
boundary- layer/ shock  interaction,  and  high  heat-transfer 
rates.  This  study  concerns  the  inviscid  flow  about  blunt 
bodies  in  the  absence  of  chemical  reactions. 

The  very  nature  of  high  velocity  and  high-temperature 
fluid  flow  makes  it  difficult  and  expensive  to  analyze  or  test 
experimentally.  Until  recently,  the  mathematical  non- 
linearities  associated  with  detached  and  curved  shocks,  as 
well  as  high-temperature  gas  interactions  and  reactions,  have 
prevented  a  numerical  approach  to  studying  hypersonics. 
However,  with  advances  in  Computational  Fluid  Dynamics  (CFD) , 


numerical  solution  techniques  have  become  possible,  and  in 
fact  preferable. 

The  neglect  of  high-temperature  gas  effects  still  limits 
the  applicability  of  most  numerical  algorithms  to  hypersonic 
problems.  To  illustrate,  consider  a  reentry  vehicle  traveling 
at  6  km/s.  It  will  experience  temperatures  near  the  stagnation 
point  of  about  6,000  °K  for  a  reacting  gas,  but  would  endure 
an  astonishing  (and  unrealistically  high)  15,000  °K  in  a 
calorically  perfect  gas  (1:19) .  This  simplistic  example  shows 
the  importance  of  accurately  modeling  high-temperature  flows. 

In  the  rest  of  this  chapter,  the  current  status  of  both 
CFD  and  hypersonic  aerodynamics  will  be  highlighted,  as  they 
apply  to  this  study. 

A.  Historical  Background 

Numerical  Analysis  of  Blunt-Bodv  Flows 

The  mixture  of  subsonic  and  supersonic  regions  in  blunt- 
body  flows  make  consistent  analysis  difficult  to  achieve.  In 
fact,  until  1966,  no  practical  blunt-body  solution  existed  to 
provide  valid  initial  conditions  for  a  method  of 
characteristics  solution  in  the  supersonic  region.  This 
changed  in  1966  when  Moretti  and  Abbett  published  the  first 
time-marching,  finite-difference  solution  of  the  governing 
unsteady  Euler  equations.  Since  the  unsteady  Euler  equations 
are  hyperbolic  in  both  subsonic  and  supersonic  regions,  they 


are  mathematically  consistent.  Although  Moretti  and  Abbett 
employed  the  Lax-Wendroff  finite-difference  technique,  the 
philosophy  remained  the  same  when  MacCormack's  explicit, 
predictor-corrector  finite-difference  method  was  later 
applied.  MacCormack's  method  has  been  widely  used  throughout 
the  1970 's  and  1980 's  (1:169).  Unfortunately,  this  method 
requires  several  points  to  capture  shocks,  imposes  strict 
limits  on  time  step  size,  and  is  not  very  robust.  The 
restriction  on  time  step  size  is  particularly  limiting,  since 
it  delays  convergence  and  induces  a  factor  of  unpredictability 
into  the  solution  process.  A  new  class  of  schemes  has 
recently  arisen  to  address  these  issues,  namely  Total 
Variation  Diminishing  (TVD)  schemes. 


Development  of  TVD/Roe  Schemes 
Many  different  schemes  and  algorithms  have  evolved  to 
solve  the  Euler  equations  (16,44,45,46).  Most  of  these 
numerical  schemes  arose  from  examinations  of  the  Riemann 
problem.  The  Riemann  problem  is  governed  by  the  1-D  inviscid 
and  frozen  flow  equation; 


dU  ^  dF 
dt  dx 


0 


(1) 


®  t  =  0 


if  X  <  0 
Ug  if  X  >  0 


where  and  Ug  are  the  left  and  right  interface  states.  By 


replacing  the  initial  data  with  an  approximate  distribution, 
which  is  constant  and  uniform  over  each  discrete  interval,  we 
can  solve  the  Riemann  problem  for  each  interval  and  get  an 
"exact"  solution  to  an  "approximate"  problem  (36:344).  This 
was  first  proposed  by  Godunov  and  is  appropriately  named 
Godunov's  Method.  As  one  might  expect,  this  approach  places 
the  discontinuities  at  cell  interfaces  (35:358).  Most  of  the 
upstream  differencing  schemes  make  use  of  the  solution  of  the 
Riemann  problem  as  a  building  block.  These  Riemann  solvers 
divide  the  flux  difference  between  neighboring  states  into 
component  parts  associated  with  each  wave  field  (11)  and  has 
yielded  great  success  in  flow  computation. 

From  a  numerical  point  of  view,  however,  it  seems 
wasteful  to  solve  exactly  the  Riemann  problem  at  every 
interface.  Therefore,  some  approximate  solutions  or  solvers 
are  preferred.  This  analysis  uses  Roe's  method  of  "... 
exploiting  the  fact  that  the  Riemann  solution  for  any  set  of 
linear  conservation  laws  is  easily  computed"  (36:346).  Here 
we  replace  F  in  Equation  (1)  with  some  linear  function  of  U  : 


dU 

dt 


dU 

dx 


=  0 


where  A  is  the  Jacobian  matrix 


dF 

du' 


(2) 


This  approach  is  much  more  efficient.  In  addition  to  the 
complexity  associated  with  solving  the  nonlinear  Euler 


equations,  it  is  difficult  to  accurately  resolve  the  shocks 
present  at  hypersonic  speeds.  For  problems  with  strong  shock 
waves  or  contact  discontinuities.  Total  Variation  Diminishing 
(TVD)  schemes  provide  high-resolution,  shock  capturing  without 
loss  of  smoothness.  In  other  words,  a  shock  is  captured  in 
only  a  few  discrete  nodes  without  sacrificing  a  smooth 
solution  near  the  discontinuity.  More  formally,  these  TVD 
schemes  comprise  a  class  of  conservative,  nonoscillatory , 
shock  capturing  schemes  (47:6).  Additionally,  when  applied  to 
a  scalar,  hyperbolic,  conservation  law,  as  (1)  or  (2),  the 
total  variation  of  a  solution,  TV{U)  ,  given  by: 


TV(U) 


(3) 


cannot  increase  with  time  (20:2).  The  first-order  Roe  scheme, 
which  was  used  in  this  analysis,  satisfies  the  TVD  property 
for  a  scalar,  hyperbolic  equation. 

Roe,  who  originally  developed  his  algorithm  in  1981,  used 
information  from  two  neighboring  states  to  compute  quantities 
at  some  intermediate  state  (35:357-372),  much  like  simple 
straight  averaging.  Unlike  straight  averaging,  however.  Roe¬ 
averaging  is  mass  weighted.  Vinokur  and  Liu  recently  extended 
Roe's  work  for  equilibrium  and  nonequilibrium  flows  (24,42). 
Chapter  2  provides  some  details  about  this  averaging. 
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B.  Hiah-Teg.perature  Effects 

Although  many  robust  techniques  use  TVD  schemes  and 
approximate  Riemann  solvers,  they  still  assume  an  ideal  gas 
and  neglect  high-temperature  effects  such  as  molecular 
vibration,  dissociation,  and  ionization.  It  is  well  known 
that  high-temperature  effects  play  an  important  role  at 
hypersonic  speeds.  Particularly,  the  shock  layer  decreases  in 
thickness  as  temperature  behind  the  shock  increases.  This  is 
caused  by  processes,  other  than  the  normal  translational  and 
rotational  degrees  of  freedom  in  molecules,  absorbing  energy 
behind  the  shock.  The  first  of  these  processes  to  develop,  as 
temperature  increases,  is  molecular  vibration.  When 
vibrational  excitation  absorbs  energy,  the  temperature 
decreases.  The  decrease  in  temperature  accompanies  a  rise  in 
density,  a  decrease  in  Y/  and  ultimately  results  in  a  thinner 
shock  layer  (31:2).  To  illustrate,  the  following  relation 
between  density  and  shock  standoff  for  an  axisymmetric  blunt- 
body  is  provided  from  Hayes  and  Probstein  (12:160): 


=  e  R. 


1  - 


8e 

3 


.|e 


+  0(6^) 


(4) 


Where  e  =  p^/pa-  Clearly,  as  P2  increases,  Aq  decreases.  Now 
by  applying  the  hypersonic  limit  to  the  normal  shock  relation 
for  density,  e  becomes: 


lim  e  =  -Xzi. 

Y+i 


6 


This  simple  expression  depends  entirely  on  the  value 
of  Y  within  the  shock-layer,  and  therefore  the  thermal  state 
of  the  flow.  The  trend  in  shock  standoff  distance  with 
decreasing  y  can  be  seen  in  Figure  1. 


Figure  1  Theoretical  Shock  Standoff  Distance 


Much  effort  recently  has  been  directed  toward  modeling 
thermal  nonequilibrium  phenomena.  In  fact,  current  schemes 
employ  various  vibration  modeling  techniques  (5,9,27,31,39). 
However,  this  modeling  is  most  often  done  in  conjunction  with 
chemical  reactions.  Since  vibrational  relaxation  precedes 
dissociation,  and  the  dissociation  rate  depends  significantly 


on  the  vibrational  state,  an  accurate  analysis  of  molecular 
vibration  is  required  in  isolation  from  the  chemistry. 

C.  Outline  of  Study 

This  study  has  three  phases.  In  the  first  phase,  the 
effects  of  high-temperature  thermodynamic  states 
(equilibrium/nonequilibrium)  are  investigated.  Thermal 
equilibrium  is  analyzed  by  using  variable  heat  capacities. 
The  equilibrium  analysis  requires  no  change  to,  or  expansion 
of,  the  equations  of  motion.  The  nonequilibrium  analysis  is 
accomplished  by  expanding  the  governing  equations  to  include 
a  vibrational  energy  equation.  The  rate  of  energy  exchange 
between  translational  and  vibrational  energy  modes  will  be 
examined  using  two  different  relaxation  rate  models.  Chapter 
3  details  the  theory  associated  with  each  model.  Computed 
results  are  compared  to  experimental  results  for  a  shock  tube. 

During  the  second  phase,  a  Roe  scheme  was  modified  for 
axisymmetric  flow  and  a  general  coordinate  system.  The  "Roe- 
averaged”  state  was  derived  by  Roe  for  a  perfect  gas  (35:357- 
372)  and  is  outlined  in  Chapter  2.  This  phase  established  a 
perfect  gas  baseline  for  comparison  throughout  the  study. 
Solutions  are  computed  over  the  forward  portion  of  a  sphere 
for  Mach  numbers  ranging  from  5  to  8.  Theoretical  shock 
standoff  distances,  sonic  line  locations  and  pressure 
distributions  are  used  to  validate  the  algorithm. 


For  the  third  and  final  phase,  contrasts  of  accuracy  and 
computational  complexity  highlight  the  advantages  of  each 
model.  This  gives  future  investigators  information  about  the 
tradeoffs  associated  with  vibrational  excitation  modeling. 

D.  Findings 

The  first-order  Roe  scheme  accurately  predicted  inviscid 
flow  about  a  sphere.  Additionally,  large  Courant  numbers 
could  be  used  during  the  analysis.  Converged  thermal 
nonequilibrium  solutions  were  obtained  in  less  than  1000 
iterations. 

Vibrational  excitation  caused  shock  standoff  distances  to 
decrease  by  about  10-20  percent.  Subsequently,  temperatures 
within  the  shock  layer  and  on  the  body  surface  differed 
substantially  from  perfect  gas  calculations.  Only  the 
nonequilibrium  analysis  predicted  a  concentrated  region  of 
high  temperatures  near  the  shock. 

Finally,  this  study  revealed  that  using  a  thermal 
nonequilibrium  model,  with  a  Millikan  and  White  relaxation 
rate,  provided  the  best  combination  of  accuracy  and 
efficiency. 


II.  Mathematical  Model 


Two  different  geometries  were  analyzed  in  this  study,  the 
1-D  shock  tube  and  a  simple  axisymmetric  blunt-body.  First, 
the  shock  tube  allowed  comparisons  between  computed  solutions 
and  theory  to  ascertain  the  performance  of  the  algorithm  in 
one  dimension  and  the  accuracy  of  each  molecular  vibration 
model.  Also,  the  relaxation  phenomena  could  be  analyzed 
without  the  complicating  effects  of  grid  transformations  and 
geometry.  The  blunt-body  problem  provided  a  more 

comprehensive  check  on  the  accuracy  in  computing  flows  around 
more  physically  realistic  geometries.  As  pointed  out  by 
Reference  21,  "at  hypersonic  flight  speeds  all  bodies  must  be 
blunt-nosed  to  some  extent  to  reduce  the  heat  transfer  rate  to 
manageable  proportions." 

A.  Assumptions 

To  fully  isolate  the  kinetic  behavior,  the  following 
assumptions  were  made  throughout  the  investigation: 

1.  Inviscid,  nonconducting  flow 

2.  All  energy  modes  decouple  completely 

3.  No  electronic  effects 

4.  No  dissociation  or  ionization 

The  decoupling  of  energy  modes  follows  from  an  analysis  of  the 
Schrodinger  equation.  This  equation  is  a  time-dependent 
equation  which  governs  two  interacting  particles. 
Translational  energy  decouples  mathematically  by  a  simple 


separation  of  variables  in  the  Schrodinger  equation  (15)  . 
Since  the  internucleaic  spacing  closely  matches  the  classical 
equilibrium  spacing  (15) ,  a  further  decoupling  of  rotation  and 
vibration  results.  A  more  detailed  discussion  of  the  quantum 
mechanics  of  diatomic  molecules  can  be  found  in  Reference  15. 

Although  viscous  drag  and  heat  transfer  cannot  be 
determined,  a  comparison  of  density  and  temperature  within  the 
shock  layer  can  be  made.  Ultimately  these  quantities  will 
influence  the  size  of  gradients  within  a  boundary-layer. 

Some  further  assumptions  for  equilibrium  air  were  made: 

1.  Vibrational  mode  in  equilibrium  with  translation  and 
rotation,  but  not  necessarily  fully  excited. 

2.  The  flow  can  be  described  by  a  single  density  and 
temperature . 

3 .  A  harmonic  oscillator  approximation  for  the  molecular 
potential  energy. 

4.  A  characteristic  vibrational  temperature,  0^,,  of 
3154.8  °K. 

The  term  "fully  excited"  means  that  the  specific  internal 
energy  increases  linearly  with  T  (41:133).  This  concept  is 
used  by  anyone  assuming  a  calorically  perfect  gas.  A  single 
density  and  temperature  assumption  simplifies  the  governing 
equations,  and  reduces  the  rank  of  the  system  to  4.  Moreover 
a  single  harmonic  oscillator  can  be  used  in  the  vibrational 
excitation  analysis,  A  harmonic  oscillator  potential  is  a 
classical  linear  spring  potential  where  the  spring  constant 
relates  to  the  vibrational  frequencies.  In  practice,  these 
frequencies  are  determined  from  spectroscopic  measurements. 
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Harmonic  vibration  reasonably  predicts  the  vibrational  energy 
states  associated  with  quantum  mechanics  (15) .  Again  the 
reader  is  referred  to  Reference  15  for  a  more  thorough 
discussion  of  quantum  and  statistical  mechanics. 

No  characteristic  vibrational  temperature  explicitly 
exists  for  air.  Therefore,  those  associated  with  oxygen  and 
nitrogen  were  used  to  derive  a  mole-based  6„  .  These  two 

''•ix 

constituents  were  selected  since  they  comprise  the  primary 
diatomic  species  in  air. 

A  harmonic  oscillator  potential,  and  the  same 
characteristic  vibrational  temperature,  were  also  used  in  the 
nonequilibrium  analysis.  Here,  however,  a  different  equation 
governs  the  energy  associated  with  molecular  vibration.  This 
equation  is  provided  in  the  next  section. 

B.  Governing  Equations 

1-D  Shock  Tube 

The  flow  in  an  ideal  shock  tube  is  caused  by  the 
instantaneous  rupture  of  a  dividing  wall,  or  diaphragm,  that 
separates  two  fluids  with  different  static  pressures.  A 
schematic  of  the  flow  field  is  as  follows: 


Figure  2  Shock  Tube  Schematic 


For  this  analysis,  the  region  between  the  shock  wave  and 
contact  surface  provide  the  proper  physical  environment  for 
nonequilibrium  study.  Additionally,  the  shock  tube  is  a 
Riemann  problem  which  has  an  exact  solution.  The  normal  shock 
relations  combine  with  Riemann  invariants  to  produce  the 
following  expression  involving  the  pressure  jump  across  the 
shock: 


E±  =  El. 

Pi  Pi 


(Y4-I)  p 

(4y^2Yi(Yi+1)P)"^" 


yZy./CYj-l) 

> 


(5) 


where,  P  =  —  -  1,  represents  the  speed  of  sound  in  front 

Pi 

of  the  moving  shock,  and  the  speed  of  sound  in  front  of  the 
expansion  fan.  The  remainder  of  the  exact  solution  is 


outlined  in  Reference  8. 


For  the  shock  tube,  the  continuity,  momentum,  and  energy 
equations  can  be  written  in  conservation  form  as: 


dU  ^  dF 
dt  dx 


0 


(6) 


where 
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Sphere 

The  equations  of  motion  for  inviscid,  non-conducting, 
axisymmetric  flow  can  be  written  in  conservation  form  as: 


where 


dU  dF  dG  _ 

dt  dx  dy 
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Appendix  A  transforms  and  develops  these  equations  for  a 
general  coordinate  system  in  the  manner  described  in  Reference 
38.  When  discretized,  the  equations  of  motion  are  evaluated 


at  cell  centers.  Additionally,  at  cell  interfaces,  the 
metrics  were  averaged  using  a  finite  volume  approach  (46:4). 

Eguilibrium  analysis  does  not  require  any  modification  of 
the  equations  of  motion.  For  the  nonequilibrium  analysis, 
however,  the  equations  of  motion  must  be  expanded  to  include 
a  separate  equation  for  the  vibrational  energy  of  each 
species.  Now  U  in  Equation  (5)  becomes: 
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^vib,y 
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The  vibrational  energy  equation  written  in  conservation  form 
is  (5:2)  . 


at 


dx 


djvE^ity'i 

dy 
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(9) 


where , 


Qj._y  is  the  Landau-Teller  form  for  the  rate  of  energy  exchange 
between  translational  and  vibrational  energy  modes,  and  x  is 
the  relaxation  time  (19) .  The  Landau-Teller  rate  equation 
assumes  that  the  vibrational  relaxation  rate  varies  linearly 
with  the  difference  in  vibrational  energies  (32:235) .  This  is 
the  most  widely  used  form  for  vibrational  relaxation  modeling 


(4,5,9,16,27,33,40,43).  TWO  different  relaxation  rate  models 
were  used  to  predict  x .  Chapter  3  describes  the  theory  behind 
each  model. 

The  linearized  flux  Jacobians  used  in  the  equilibrium  and 
nonequilibrium  analysis  were  obtained  with  Mathematica,  which 
is  a  system  for  doing  symbolic,  algebraic  and  graphic 
mathematical  computations.  Appendix  B  contains  a  summary  of 
the  procedure  and  resulting  eigenvalues  and  eigenvectors. 


C.  The  Roe-averaaed  State 

When  evaluating  the  dissipation  term,  some  form  of 
averaging  is  required  at  cell  interfaces.  The  present  scheme 
uses  the  "Roe-average".  This  Roe-averaged  state  was  derived 
by  Roe  for  a  perfect  gas,  and  generalized  for  both  equilibrium 
and  nonequilibrium  flow  by  Vinokur  (24,35,42).  The  reader  is 
referred  to  these  references  for  a  complete  derivation  of  the 
averaging.  To  summarize.  Roe-averaging  consists  of, 
^i*i/2  ^i*i/2  -  f  makes  use  of  the 

ratio  Pj.i/Pi  such  that: 


.1/2 


1/2. 
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Pi  Ui  +  PiViUi^i 
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-1/2 
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(10) 


where,  D  =  (pj^^/p^)  .  Since  pressure  is  directly  related  to 
enthalpy  for  a  calorically  perfect  gas,  typically  only 
velocity,  enthalpy  and  the  speed  of  sound  need  to  be  averaged. 
However,  for  the  equilibrium  and  nonequilibrium  analysis 


pressure  no  longer  directly  relates  to  enthalpy  and  the 
partial  derivatives  of  pressure,  defined  here  as  k  and  %/ 
also  need  to  be  averaged.  A  detailed  discussion  of  these  two 
parameters  can  be  found  in  the  next  chapter,  but  the  averaging 
for  these  quantities  is  quite  involved  and  will  not  be 
included  into  this  document.  References  24,42  amply  cover  the 
derivations  in  full  detail. 

D.  Initial  and  Boundary  Conditions 

The  shock  tube  problem  only  requires  that  the  initial 
solution  have  a  high  pressure  side  and  low  pressure  side.  It 
does  not  require  any  boundary  conditions  unless  reflected 
shocks  and  waves  are  studied. 

For  the  blunt-body  problem,  an  initial  condition  was 
derived  from  an  empirical  shock  location  and  shape,  oblique 
shock  relations,  modified  Newtonian  theory,  the  equation  of 
state  and  a  constant  enthalpy  condition  (1:190).  The  shock 
standoff  distance  and  shock  shape  can  be  very  closely 
approximated  by  the  method  in  Reference  2.  Oblique  shock 
relations  then  determine  values  at  the  shock.  Along  the  body, 
the  initial  condition  assumes  a  Newtonian  pressure 
distribution  and  linear  velocity  distribution.  Throughout  the 
field,  velocities  and  pressures  are  interpolated  between  the 
body  and  shock.  Finally,  the  adiabatic,  constant  enthalpy 
relation  and  equation  of  state  determine  temperature  and 
density,  respectively. 


Boundary  conditions  which  were  simple,  but  would  not 
degrade  the  accuracy  of  the  solution  were  chosen  for  this 
study.  For  inviscid  flow,  tangency  of  velocity  needs  to  be 


satisfied  at  the  body  surface.  Reflection  on  density  and 
total  energy  complete  the  surface  boundary  condition.  Next, 
a  supersonic  inflow  condition,  and  "no-change”  outflow 


condition  were  imposed. 
dU 


Finally,  at  the  symmetry 


line  V  and  -*  0  . 
dr 


Several  other  boundary  conditions  were 
tried  while  trying  to  improve  the  solution: 


BODY 

1.  Tangency  and 

Normal  Momentum  (34:14) 


SYMMETRY  LINE 

1.  Zeroing  the  fluxes 

2.  Limiting  form  of  the 
Equations  of  Motion(39:4) 


Chapter  4  summarizes  results  from  the  different  conditions. 
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Ill,  vibrational  Excitation  Theory 


Vibrational  excitation  can  be  modeled  in  three  ways.  We 
can  ignore  it,  as  with  the  perfect  gas  assumptions.  We  can 
also  place  it  in  equilibrium  with  the  other  energy  modes,  or 
we  can  treat  it  as  a  process  that  occurs  at  a  finite  rate.  As 
described  by  Lambert  (18:13): 


In  a  sudden  compression  the  whole  energy  goes 
initially  to  increase  and  this  is  followed 
by  relaxation  to  T-R  and  T-V  energy  transfer 
until  equilibrium  between  the  three  modes  is 
re-established.  Such  rapid  compressions  are 
produced  by  ultrasonic  sound  or  by  the  passage 
o f  shock- f rents . 


A.  Equilibrium  Vibrational  Energy 

The  equilibrium  model  is  based  solely  on  statistical  and 
quantum  mechanics,  which  predicts  that  a  partition  function 
can  be  related  to  specific  internal  energy  by  the  relation 
(41:128) : 


'vib 


=  RT^ 


a  (In 


dT 


(11) 


where  is  the  Partition  Function,  and 


^vib 


_ 1 _ 

1  -  exp(-0yr) 


Recall  that  0^  is  the  characteristic  vibrational  temperature 
and  is  constant.  Additionally,  since  internal  energy  defines 


temperature,  the  above  relation  for  the  vibrational  energy, 
combined  with  the  relations  for  fully  excited  translational 
and  rotational  energy,  provide  enough  information  to  compute 


temperature . 


=  ^RT 


^rot  =  Rl" 


Only  for  temperatures  which  exceed  those  where  dissociation 
begins  will  the  vibrational  energy  mode  reach  a  fully-excited 
state. 

As  pointed  out  earlier,  the  equilibrium  model  required  no 
expansion  or  modification  of  the  governing  equations.  The 
increase  in  complexity  finds  its  way  into  the  problem  through 
the  linearization  of  the  flux  terms.  For  these  terms,  the 
partial  derivatives  of  pressure,  with  respect  to  energy  and 
density,  were  rederived.  Utilizing  the  chain  rule: 


-  de,-[dTj[dsJ-^^de^ 


where  e  =  pe.  Now  from  (11-12), 


5  ^vR 

e  =  —RT  +  - 7  , 

2  expceyr) 
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expoyr)  pR 

±PR  +  - Uil - 

2  (expe^/r  -  1)2 


exp(eyT) 

J.  + _ _ _  ^ 

2  (expeyr  -  1)2 


Using  the  definition  of  pressure  from  Reference  42,  leads  to 
a  convenient  form  for  the  second  partial  derivative  of 
pressure: 

P  =  Xp  +  K§ 


dp 

where  X  = 

dp  g 


Therefore, 


X  =  RT  -  ^ 


It  is  very  interesting  to  note  that  (16)  reduces  to 
2 

X  =  —  or  (y-l)  for  low  temperatures.  Then  with  S  =  pC^T 
5  ^ 

and  K  =  Y~l/  X  from  (17)  vanishes.  These  are  exactly  the 
perfect  gas  relations  where  p  =  (Y-l)pe. 
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B.  Nonecniillbrium  Relaxation 


Equilibrium  for  translation  is  attained  in  only  a  few 
collisions  since  there  is  free  interchange  of  translational 
energy  at  every  collision.  Likewise,  for  temperatures  in 
excess  of  a  few  hundred  degrees,  rotational  equilibrium  is 
attained  in  less  than  ten  collisions  (18:25).  Molecular 
vibration,  however,  requires  significantly  more  collisions  and 
time  to  reattain  equilibrium.  The  present  nonequilibrium 
analysis  accounts  for  this  energy  relaxation  with  a  two- 
temperature  model  (30:488).  One  temperature,  T,  specifies 
the  translational/rotational  energy,  and  the  other 
temperature,  defines  the  nonequilibrated  vibrational 
energy.  This  analysis  used  the  Landau-Teller  form  for  the 
energy  exchange  between  modes  (19:34).  This  appears  as  a 
source  term  in  the  vibrational  energy  equation  (see  Chapter 
2),  and  is  a  function  of  the  relaxation  time,  x.  To 
approximate  the  relaxation  time,  both  the  Millikan  and  White 
model  and  the  Park  model  were  tested.  These  models  are  both 
emperical,  since  a  theoretical  model  has  not  been  developed  to 
date. 


Millikan  and  White  Model 

Millikan  and  White  developed  an  empirical  model  to 
compute  relaxation  times  for  diatomic  molecules (28 : 1) .  They 
found  that  In(xp)  is  linearly  related  to  through  the 


following  relation: 


where , 


=  -^expU^^(r-i/3_|^i/4)  _^g  42]  (18) 

^  S 

=  1.16x10-3^^/2  04/3 

=  K  ^X) 


Here  Xg^  is  the  relaxation  time  for  molecule  s  colliding  with 
species  r,  and  is  the  partial  pressure  of  the  colliding 
particles  (in  atm) .  As  Mach  number  increases,  both  p  and  T 
within  the  shock  layer  increase.  This  forces  relaxation  times 
to  decrease. 


Park  Model 

Park  developed  another  model  to  improve  relaxation  rates  at 
very  high  temperatures  (31:4).  The  rate  of  change  of 
vibrational  energy  in  the  molecule  s,  by  collisions  with 
species  r,  is  given  by: 
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The  number  density  of  the  mixture  is  n.  A:  is  the  Boltzsmann 

► 

constant;  jT  defines  the  average  mass  and  denotes  the 
temperature  behind  the  shock.  c  is  the  average  molecular 
speed  and  a  is  the  limiting  cross-section.  Park  modified  the 

► 

classical  Landau-Teller  model  to  account  for  the  collision 
time,  t^.  Unfortunately,  he  did  not  document  the  development 
of  the  last  term,  other  than  to  state  that  it  brought 
"calculated  values  into  closer  agreement  with  experimentally- 
determined  values"  (31:4).  The  vibrational  energy,  is 

still  determined  using  a  harmonic  oscillator  for  the  partition 
function.  The  dependence,  however,  is  now  on  the 

nonequilibrated,  vibrational  temperature,  Ty^j^  : 

I  ^ 

f  p  s  _ (9Q\ 

exp(e,yr,,^)  -  1 

►  Although  this  model  yielded  results  closer  to  experimental 

values  in  Reference  31,  no  information  on  its  application  at 
moderate  Mach  numbers  could  be  found. 
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IV.  High-Temperature.  Shock  Tube  Results 


A.  Perfect  Gas 

The  upwind  scheme  with  Roe  averaging  was  very  robust  and 
produced  gualitatively  accurate  solutions.  Figures  3-6 
demonstrate  the  accuracy  of  the  algorithm  developed.  Here  the 
distributions  for  density,  pressure,  velocity,  and 
temperature,  are  plotted  against  the  exact,  theoretical 
solution. 


Figure  3 


Pressure  Profiles  for  Perfect  Gas  Model  Applied 
to  Shock  Tube 


Density  (kg/m^)  ^  Velocity  (m/s) 


e  4  Velocity  Profiles  for  Perfect  Gas  Model  Applied 

to  Shock  Tube 
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Figure  6 


Temperature  Profiles  for  Perfect  Gas  Mod;il 
Applied  to  Shock  Tube 


The  scheme  predicted  with  reasonable  accuracy  all  three  flow 
phenomena  (shock,  expansion,  and  contact  discontinuity) . 
Furthermore,  the  scheme  smoothly  captured  the  moving  shock  in 
only  three  or  four  points,  which  is  quite  excellent  for  the 
strength  of  the  shocks  involved.  First-order  accuracy, 
however,  is  apparent  from  the  errors  in  the  expansion  fan  and 
contact  discontinuity. 

An  interesting,  nonphysical  solution  appeared  in  the 
shock  tube  numerical  simulations.  Figure  7  shows  a 
discontinuity  within  the  expansion  fan  that  does  not  arise  in 
theory  or  experiment. 
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Figure  7  Nonphysical  Discontinuity  in  Shock  Tube 


This  erroneous  phenomenon  was  eliminated  by  satisfying  the 
"entropy  condition  for  shocks"  (10:297-305).  The  entropy 
condition  is  needed  to  pick  out  the  physically  relevant 
solution,  because  solutions  of  hyperbolic  conservation  laws 
are  not  uniquely  determined  by  their  initial  values  (10:297). 
A  small  constant  parameter,  as  proposed  by  Reference  47,  was 
incorporated  into  the  dissipation  term  and  proved  sufficient 
to  prevent  nonphysical  solutions. 


B.  Vibrational  Effects 


Figures  8-12  show  the  shock  tube  computations  using  each 
vibrational  model.  The  pressure  is  only  nominally  effected  by 
vibrational  excitation,  and  changed  by  about  two  percent  from 
theoretical  predictions  (Figure  8) . 


Figure  8  Pressure  Profile  Comparison  for  Shock  Tube 


Temperature  and  density,  however,  deviated  significantly  (10 
15  percent)  from  perfect  gas  theory  (Figures  10  and  11) . 
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e  9  Velocity  Profile  Comparison  for  Shock  Tube 
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Figure  10  Density  Profile  Comparison  for  Shock  Tube 
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Figure  11 


Temperature  Profile  Comparison  for  Shock  Tube 
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Figure  12  clearly  shows  the  nonequilibrium  relaxation 
phenomena.  The  translational/rotational  temperature  reaches 
a  peak  value  directly  behind  the  shock,  then  quickly  falls  as 
molecular  vibration  begins  to  absorb  energy.  Equilibrium 
exists  where  T  and  '^vib  are  equivalent. 


Figure  12  Vibrational  Relaxation  in  Shock  Tube 
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Density  Ratio  Calculations 


Both  qualitative  and  quantitative  comparisons  were  made 
with  experimental  density  information.  The  qualitative 
comparisons  determined  the  ability  of  each  model  to  predict 
the  flow  features  in  the  shock  tube.  The  quantitative 
comparison  provided  a  detailed  measure  of  the  accuracy. 

To  make  a  qualitative  assessment  of  the  equilibrium  and 
nonequilibrium  models,  density  interferograms  were  created 
from  the  computed  densities  at  two  different  conditions. 
These  were  then  compared  to  experimental  interferograms.  The 
first  comparison  was  for  a  shock  moving  at  a  low  Mach  number 
(M=2) .  For  this  case  minimal  molecular  vibration  should  take 
place.  The  second  comparison  was  for  a  shock  moving  at  a  much 
higher  Mach  number  (M=6) .  Here  molecular  vibration  is 
significant.  The  interferogram  is  an  interference  measurement 
which  shows  jumps  and  gradients  in  the  flow  field.  Gradients 
appear  as  curvature  in  the  interference  pattern. 

Figure  13  (a,b,c)  shows  the  interferograms  for  M=2.  Both 
the  equilibrium  and  nonequilibrium  models  match  the  experiment 
and  predict  a  nearly  constant  density  behind  the  shock. 
Figure  14  (a,b,c)  shows  the  interferograms  for  M=6.  Only  the 
nonequilibrium  model  predicts  the  relaxation  phenomenon  shown 
in  the  experimental  case. 


•  Figure  13a  Experimental  Density  Interfsrogram 

(Shock  Tube,  M=2) 


Figure  13b  Nonequilibrium  Density  Interferogram 
^  (Shock  Tube,  M=2) 


• 

Figure  13c 

Equilibrium  Density  Interferogram 
(Shock  Tube,  M=2) 
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• 

Figure  14a 


Experimental  Density  Interferogram 
(Shock  Tube,  M==6) 


Figure  14b  Nonequilibrium  Density  Interferogram 
i  (Shock  Tube,  M=6) 
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Next  more  quantitative  comparisons  were  made  to  fully 
assess  each  models  validity.  To  provide  a  common  point  of 
comparison,  the  densities  were  calculated  after  the  flow  re¬ 
equilibrates  in  the  experimental  and  nonequilibrium  cases. 
The  ratios  across  a  normal  shock,  for  both  nonequilibrium 
(relaxing)  and  equilibrium  models,  matched  the  experimental 
results  (3:75-79)  for  nondissociating  air  within  two  percent. 
Table  I  summarizes  results  for  a  shock  wave  moving  at  M=5.75. 
Equilibrium  may  match  the  experimental  results  better  due  to 
experimental  error. 


Model 

^2 

P2 

Pz 

Pi 

Pi 

Theory 

7.37 

5.21 

38.41 

Perfect  Gas 

7.36 

5.20 

38.36 

Equilibrium 

7.00 

5.80 

40.55 

Nonequi 1 ibri um 

6.89 

5.84 

40.03 

Experiment 

5.75 

Table  I  Calculations  across  a  Moving  Shock  (M  =  5.75). 


Experimental  density  ratios  for  oxygen  and  nitrogen  are 
plotted  in  Figure  15  with  the  results  from  the  nonequilibrium 
scheme  using  the  Millikan  and  White  model.  Clearly,  the 
computations  for  air  follow  the  experimental  trend. 
Additionally  the  computed  ratios  for  air  fall  between  oxygen 
and  nitrogen,  as  one  would  expect.  Recall  that  the 


I 


characteristic  vibrational  temperature,  0^,,  used  for  air 
(3155  °K  )  was  a  mole-based  average  of  0,,  for  oxygen (2270  °K  ) 


and  0^,  for  nitrogen  (3390  )  . 


Figure  15  Density  Ratios  Across  a  Moving  Shock 


Relaxation  Times 

►  The  relaxation  time  was  computed  by  observing  a 

stationary  point  in  the  shock  tube.  After  the  numerical  shock 
passed  over  this  point,  time  was  measured  until  the  flow 
•  reattained  equilibrium.  Figure  16  shows  a  comparison  of 
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relaxation  times  for  pure  nitrogen  and  pure  oxygen,  normalized 
to  1  atm.  Relaxation  times  compared  favorably  for  the  range 
of  moderate  temperatures  (Mach  numbers)  tested.  Although 
exact  agreement  was  not  obtained,  quantities  computed  with  the 
Millikan  and  White  model  fall  close  to  the  scatter  of 
experimental  data. 
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Figure  16  Oxygen  and  Nitrogen  Relaxation  Times  normalized 
to  1  atm. 


After  agreement  for  pure  species  was  confirmed,  mole- 
averaged  constants  were  employed  in  the  calculations  for  air 
as  a  single,  composite  species.  Figure  17  shows  a  comparison 


» 

of  relaxation  times  between  experiment  and  the  nonequilibrium 
scheme  for  air  with  the  Millikan  and  White  model.  As  can  be 
^  seen,  at  temperatures  below  2500  degrees  K,  the  scheme 

predicts  relaxation  rates  between  pure  oxygen  and  pure 
nitrogen;  much  closer  to  nitrogen.  This  is  consistent  with 
^  the  mole-averaged  0^  used.  At  higher  temperatures,  the 


relaxation  rates  are  nearly  identical  to  those  of  nitrogen. 


Figure  17  Relaxation  Times  for  Air,  Normalized  to  1  atm. 


Two  possible  explanations  for  the  models  inaccuracy  at 
higher  temperatures  are  a  slightly  high  characteristic 
vibrational  temperature  or  a  failure  to  account  for  the 
dissociative  behavior  of  oxygen  at  these  temperatures.  A 


characteristic  vibrational  temperature  for  air  does  not 
explicitly  exist.  The  mole-averaged  quantity,  derived  for 
this  analysis,  does  not  take  into  account  the  effectiveness  of 
different  collision  types,  like  O2-O2  and  O2-N2.  This 
significantly  impacts  relaxation  rate  prediction  as  pointed 
out  by  Reference  2.  Furthermore,  oxygen  will  begin  to 
dissociate  above  2300  degrees  K. 


Park  Model  Versus  Millikan  and  White  Model 
Only  minor  differences  between  the  two  models  could  be 
noted  for  the  temperature  range  investigated.  Table  II  gives 
a  comparison  of  results  using  the  Park  model  and  the  Millikan 
and  White  model. 


Temperature 
fdea  K) 

Park  Model 
fmicrosec^ 

Mil.  &  White  Model 
rmicrosec) 

3550 

21.4 

16.4 

3350 

27.3 

21.3 

3200 

34.1 

26.1 

3050 

40.8 

32.2 

2950 

47.4 

37.8 

2850 

53.9 

44.1 

2750 

59.0 

51.5 

2600 

69.4 

65.2 

2400 

84.3 

80.8 

2200 

186.7 

186.5 

2050 

257.8 

257.7 

1850 

577.0 

577.0 

1800 

668.2 

668.2 

Table  II  Relaxation  Times  for  Two-Temperature  Models. 


As  temperature  increased,  the  Park  model  calculated  longer 
relaxation  times  than  the  Millikan  and  White  Model.  This  is 
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consistent  with  the  design  of  this  model.  As  stated  by 
Park(32:234-235)  : 


The  modification  accounts  for  two  phenomena 
that  are  unique  to  the  high-temperature 
e  .vironment,  that  is,  a  limiting  value  of  the 
cross  section  and  the  diffusive  nature  of 
vibrational  relaxation.  The  first  phenomenon 
arises  from  the  fact  that  the  correlation 
formula  of  Millikan  and  White  implies  an 
unrealistically  large  cross  section  for 
vibrational  relaxation  at  high  temperatures. 
The  second  correction  accounts  for  the  fact 
that,  at  high  temperatures,  vibrational 
relaxation  obeys  a  diffusion  equation  with 
respect  to  the  vibrational  energy  levels, 
rather  than  a  Landau-Teller-type  rate 
equation.  As  a  result,  the  vibrational 
relaxation  rate  does  not  vary  linearly  with 
the  difference  in  vibrational  energies. 


The  diverging  trend,  as  temperature  increases,  is  visible  in 
Figure  18. _ 
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Figure  18 


Comparison  of  Relaxation  Times  for  Two 
Temperature  Models. 


No  other  computational  results  for  the  temperature  range 
investigated  could  be  found  to  confirm  these  results.  However 
Park  developed  his  model  for  dissociating  and  ionizing  flows 
where  the  mean  temperature  exceeds  10,000  °K.  These  results 
indicate  that  the  Millikan  and  White  model  is  more  accurate  at 
temperatures  below  4000  °K.  Therefore,  for  the  remainder  of 
this  study,  only  the  Millikan  and  White  Model  was  used. 


V.  Blunt-Bodv  Solutions  for  a  Perfect  Gas 


Solutions  of  the  Euler  equations  are  presented  in  this 
section  for  perfect  gas  flows  about  the  forward  portion  of  a 
sphere.  Shock  standoff  distances  agreed  with  experimental 
results,  and  sonic  line  locations  agreed  with  theory.  The 
first-order  Roe  scheme  was  very  robust  to  both  starting 
condition  and  time  step  while  calculating  blunt-body 
solutions.  Unfortunately,  the  algorithm  was  sensitive  to  the 
sphere's  geometric  singularity  and  tended  to  produce 
oscillations  and  other  nonphysical  solutions  in  the  vicinity 
of  the  stagnation  point.  Additionally,  some  sensitivity  to 
grid  refinement  was  noted. 

A.  Solution  Flow  Structure 

Qualitatively,  the  scheme  predicted  shock  standoff 
distance,  sonic  line  location,  and  surface  pressure 
distribution  quite  well.  Figure  19  shows  a  comparison  between 
experimental  and  computed  shock  standoff  distances.  Tabulated 
shock  standoff  distances  were  taken  from  References  2,17  and 
29.  Computed  values  differed  from  experiment  by  about  ten 
percent.  This  is  primarily  due  to  the  absence  of  a  boundary- 
layer  in  this  analysis.  A  boundary-layer  would  essentially 
increase  the  effective  body  radius  and  produce  larger  shock 
standoff  distances.  The  trend  with  increasing  Mach  number  is 
nicely  predicted. 


Figure  19  Computed  versus  Experimental  Shock  Standoff 
Distance  (Nonequilibrium  Model) . 


•  The  blunt -body  calculation  also  produced  accurate  sonic 
line  information.  In  axisymmetric  flow,  the  sonic  line  should 
be  obtuse  to  the  body  at  higher  Mach  numbers  (1:184) .  Figure 

•  20  shows  the  computed  sonic  line  for  a  sphere  at  a  freestream 
Mach  number  of  7.  Additionally  the  sonic  line  should  be  acute 
for  two-dimensional  flow.  Figure  21  shows  the  computed  sonic 

•  line  for  a  cylinder  also  at  Mach  7.  These  figures  indicate 
that  the  algorithm  predicts  the  proper  physical  aspects  of  the 
sonic  line  behavior. 
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Figure  20  Sonic  Line  Location  for  a  Sphere  at  M=7.0 


% 


Figure  21  Sonic  Line  Location  for  a  Cylinder  at  M=7.0. 
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Experimental  surface  pressure  coefficients  were  extracted 
from  graphical  data  in  Reference  29.  The  experimental  data  is 
plotted  in  Figure  22  with  a  modified  Newtonian  distribution. 


Figure  22  Experimental  Pressure  Coefficients  Around  the 

Forward  Portion  of  Hemisphere  Cylinder  at  M=5.7. 


Pressure  coefficients  from  the  perfect  gas  code  under 
predicted  experimental  data  by  about  five  percent.  Although 
the  Roe  scheme  does  no  better  than  modified  Newtonian,  the 
distribution  around  the  sphere  is  well  predicted  (Figure  23) , 
especially  at  the  shoulder.  The  errors  are  primarily  induced 
by  the  algorithm's  sensitivity  near  the  symmetry  line.  This 
shortcoming  is  discussed  in  the  next  section. 


Figure  23 


Pressure  Coefficients  from  Perfect  Gas 
Calculations. 


B.  Numerical  Issues 

Sensitivity  Near  Symmetry  Line 
An  oscillation  near  the  symmetry  line  developed  in  the 
solution  of  the  axisymmetric  problem.  The  two-dimensional 
problem  did  not  encounter  these  oscillations,  but  similar 
behavior  was  observed  by  Shang  and  Joshula  (39:6).  Figure  24 
shows  Mach  contours  for  flow  over  a  sphere  at  Mach  5.7;  Figure 
25  shows  similar  Mach  contours  for  flow  over  a  cylinder  also 
at  Mach  5.7.  The  solution  oscillates,  in  the  axisymmetric 
case,  for  the  first  few  degrees  from  the  stagnation  line. 


47 


Figure  24 


Mach  Contours  for  a  Sphere  at  M=5.7. 
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Figure  26  shows  the  computed  pressure  coefficients  around  a 
cylinder.  Again  no  oscillations  appeared  in  the  2-D 
computations . 


Figure  26  Surface  Pressure  Coefficients  for  a  Cylinder 
at  M=5.7. 


Studies  have  shown  that  proper  treatment  of  numerical 
boundary  conditions  has  a  major  impact  on  the  stability  of  a 
scheme  (47:105).  Unfortunately,  several  variations  in  the 
boundary  conditions  did  not  improve  or  e...:.minate  the 
oscillatory  behavior. 
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The  boundary  conditions  attempted  include: 

1.  Tangency  and  reflection  at  the  body 

2 .  Tangency  and  normal  momentum  at  the  body 

3.  Zeroing  the  fluxes  at  the  symmetry  line 

4.  Vanishing  normal  derivatives  of  conserved  variables 
at  the  symmetry  line 

5.  Utilizing  the  limiting  form  for  the  equations  of 
motion  at  the  symmetry  line. 

Several  other  techniques  fared  no  better  at  improving  the 
solution.  These  included: 

1.  Artificial,  explicit  dissipation  near  the  symmetry 
line 

2 .  Incorporating  an  entropy  parameter  to  satisfy  the 
"entropy  condition" 

Recent  work  with  second-order  Total  Variation  Diminishing 
(TVD)  schemes  has  shown  that  limiters  and  more  sophisticated 
dissipation  functions  eliminate  the  oscillatory  behavior 
(47:106).  It  is  important  to  note  that  the  oscillations  were 
confined  to  the  first  five  azimuthal  nodes,  and  not  to  the 
first  few  degrees.  As  the  grid  was  refined,  the  region 
effected  by  the  spurious  behavior  shrank.  Therefore,  with  a 
refined  grid,  the  integrated  affect  would  be  small. 

Nonphysical  Solutions 

At  the  highest  Mach  numbers  tested,  a  nonphysical 
recirculation  developed  in  the  stagnation  region  at  the  shock. 
This  phenomenon  is  shown  in  Figure  27.  As  pointed  out  by 
Reference  39,  two  of  the  four  eigenvalues  of  the  system  of 
equations  are  zero,  therefore  initial  error  will  accumulate  in 
the  iterative  process.  This  reference  characterized  the 


result  as  an  "accentuated  bulge  in  the  shock  wave  near  the 
axis  of  symmetry"  (39:6).  They  propose  that  updating  the 
initial  condition  may  eliminate  the  error.  Attempts  to  update 
conditions  within  the  stagnation  region  did  not  yield 
promising  results  in  this  study.  In  theory,  the  error  might 
be  propagated  away  from  this  region  by  using  mesh  sequencing 
or  limiters  in  the  dissipation  formulation  (34:6-12). 


Figure  27 


Nonphysical  Recirculation  Phenomenon  in 
Stagnation  Region. 


Grid  Sensitivity 


The  recirculation  phenomenon  and  the  grid  were  highly 
coupled.  As  the  grid  was  refined,  recirculation  occurred  at 
lower  Mach  numbers.  Several  permutations  in  the  grid  proved 
unfruitful,  including  grids  which  stride  the  symmetry  line  and 
grids  which  fall  on  the  symmetry  line.  By  maintaining  a 
sufficiently  large  azimuthal  grid  spacing  near  the  symmetry 
line,  and  a  sufficiently  large  Courant  number,  the  phenomenon 
could  be  avoided.  Cell  aspect  ratios  of  about  three  to  one 
worked  best.  Courant  numbers  between  0.7  and  0.9  provided  the 
necessary  dissipation  near  the  shock. 

While  maintaining  the  optimal  cell  aspect  ratios,  the 
grid  was  refined  until  the  shock  standoff  distance  converged. 
Convergence  was  determined  from  the  percentage  change  in  shock 
standoff,  and  satisfied  when  this  change  fell  below  one 
percent.  A  plot  of  shock  standoff  distance  versus  radial  node 
spacing  is  provided  in  Figure  28.  A  coarser  grid  seems  to 
underpredict  the  shock  layer  thickness. 


standoff  Distance  (x/R) 


53 


VI.  Eauilibriuin/Nonecati librium  Blunt-Bodv  Solutions 


After  validating  the  scheme  for  the  perfect  gas  model 
with  theory  and  experimental  data,  several  comparisons  to  this 
baseline  approach  proved  illuminating.  First,  molecular 
vibration  significantly  effected  shock  standoff  distance. 
Second,  an  interesting  low  density,  high-temperature  zone 
developed  immediately  behind  the  shock.  Finally,  the 
equilibrium  model  proved  far  less  robust  than  hoped. 

A.  Shock  Standoff  Distances 

For  both  the  equilibrium  and  nonequilibrium  cases, 
density  increased  and  shock  standoff  distance  decreased. 
Density  contours  for  perfect  gas  and  equilibrium  airflows  over 
the  sphere  are  presented  in  Figure  29.  A  noticeable  increase 
in  density  and  decrease  in  shock  standoff  distance  can  be 
observed  for  the  equilibrium  case.  Figure  30  shows  a  similar 
temperature  contour  comparison  between  equilibrium  and 
nonequilibrium  air.  As  expected,  the  equilibrium  model 
underpredicts  the  shock  standoff  distance.  The  vibrational 
relaxation  phenomenon  is  clearly  visible  in  Figure  30.  An 
isolated  region  of  high  temperatures  forms  directly  downstream 
of  the  shock.  As  the  vibrational  states  of  the  molecules 
bcv^'orne  excited  and  absorb  energy,  the  temperature  must 
decrease.  The  shape  of  computed  temperature  contours,  for  the 
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Figure  29  Density  Contours  for  a  Sphere  at  M= 


nuinber  increases,  molecular  excitation  increases,  and  y 
decreases.  The  nonequilibrium  calculations  produced  changes 
consistent  with  this  theory  for  both  planar  and  axisymmetric 
flow.  Figures  31  and  32  show  Mach  contours  around  a  cylinder 
(planar  flow)  for  freestream  Mach  nvimbers  of  5  and  7 
respectively.  Note  the  considerably  larger  decrease  in  shock 
layer  thickness  for  the  higher  Mach  number. 


Figure  31 


Mach  Contours  for  a  Cylinder  at  M=5.0 


Figure  32 


Mach  Contours  for  a  Cylinder  at  Mach  7.0. 


The  axisymmetric  analysis  produced  similar  results.  Here,  for 
a  freestream  Mach  number  of  6,  Aq  decreased  by  about  7 
percent.  At  a  freestream  Mach  number  of  8,  Aq  decreased  by 
over  19  percent.  If  gamma  had  decreased  the  same  amount  for 
both  Mach  numbers,  the  decrease  in  A,,  would  also  have  been 
the  same.  Table  III  summarizes  shock  standoff  results  for  a 
Mach  number  of  5.7. 
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Table  III  Shock  Standoff  Distances  for  Hemisphere  Cylinder 
at  M=5.7. 


As  described  earlier,  the  inviscid  scheme  underpredicts 
the  shock  layer  thickness  due  to  the  absence  of  a  boundary- 
layer.  Note,  however,  shock  location  followed  the  proper 
trend  as  the  various  models  were  incorporated.  Two  physical 
aspects  of  fluid  motion  support  this  trend.  First  molecular 
vibration  absorbs  energy;  therefore,  temperature  must 
decrease.  Since  the  pressure  is  only  weakly  coupled  to  the 
translational  energy  mode,  this  molecular  vibration  will 
primarily  require  the  densities  to  be  lower.  Second, 
conservation  of  mass  must  be  maintained.  Lower  densities 
require  more  volume  to  satisfy  consistent  mass  flow 
constraints.  Both  of  these  principles  effect  both  thermal 
equilibrium  and  thermal  nonequilibrium,  but  the  nonequilibrium 
case  requires  a  finite  time  and  distance  to  accomplish  the 
energy  transfer,  and  will  thus  have  thicker  shock  layers. 
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B.  Relaxation  in  Stagnation  Region 

Nonequilibriuiti  calculations  produced  striking  temperature 
and  density  gradients  in  the  stagnation  region,  with  an 
isolated  "hot  spot"  near  the  shock.  Figure  33  shows  the 
region  of  highest  temperature  located  next  to  the  shock.  This 
contrasts  with  perfect  gas  and  equilibrium  calculations  where 
the  stagnation  point  temperature  is  highest. 


Figure  34  shows  the  translational/rotational  temperature 
plotted  against  the  vibrational  temperature  along  the 
stagnation  line.  The  flow  never  quite  reattains  equilibrium. 


Figure  34  Temperatures  Along  the  Stagnation  Line. 


The  validity  of  the  two-temperature  model  is  questionable 
at  moderate  Mach  numbers,  because  the  relaxation  rates 
discussed  in  Chapter  4  differed  significantly  for  comparable 
temperature  ranges.  Figure  35  shows  the  stagnation  line 
vibrational  temperatures  for  pure  oxygen  and  pure  nitrogen. 
The  differences  throughout  the  shock  layer  are  quite  large. 
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Figure  35  Temperature  Comparison  within  the  Shock  Layer 
(Oxygen  -vs-  Nitrogen) . 


Surface  temperature  distributions,  on  the  other  hand,  did 
not  differ  significantly.  In  fact,  the  vibrational  modeling 
had  little  impact  on  the  trend  in  temperature  distribution 
around  the  sphere.  Figure  36  shows  the  translational 
temperature  around  the  sphere  for  all  three  gas  models: 
perfect  gas,  equilibrium,  and  nonequilibrium.  Discrepancies 
between  the  magnitudes  are  apparent,  but  clearly  the  geometry 
dominates  the  trend. 
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Figure  36  Surface  Temperature  Distributions  for  M=5.7. 


C.  Molecular  Vibration  Model  Tradeoffs 

The  equilibrium  model  was  nearly  as  complex  as  both  two- 
temperature  relaxing  models.  The  equilibrium  model  required 
more  CPU  time  for  two  reasons.  First,  a  Newton  iteration 
procedure  was  employed  to  determine  the  temperature  at  each 
node  point.  Second,  several  more  parameters  had  to  be  "Roe- 
averaged”  to  determine  the  eigenvalues  and  eigenvectors  of  the 
linearized  flux  Jacobians.  Computation  times  for  the  code  on 


a  Cray-YMP  are  summarized  in  Table  IV.  No  effort  was  made  to 
fully  vectorize  the  code,  so  some  increase  in  efficiency  could 
be  realized. 


Model  CPU  Time  Recmired 

Perfect  Gas  1.75  x  10‘*  sec/node  pt/time  step 

Equilibrium  9.00  x  lO*"^  sec/node  pt/time  step 

Nonequilibrium  12.30  x  10'^  sec/node  pt/time  step 

Table  IV  CPU  Requirements  for  each  Vibrational  Model 

Equilibrium  and  nonequilibrium  vibrational  modeling 
produced  similar  surface  temperature  distributions.  At  the 
lowest  Mach  numbers  tested,  the  differences  were  less  than  15 
percent.  At  Mach  8,  the  difference  shrunk  to  only  3  percent. 

Only  nonequilibrium  modeling  properly  captured  the  radial 
temperature  gradients  within  the  shock  layer.  This 
significantly  effected  proper  resolution  of  temperatures 
within  the  stagnation  region,  especially  at  or  near  the  shock. 


VII.  Conclusions  and  Recommendations 


A  first-order  Roe  scheme  was  developed  for  inviscid 
axisymmetric  flows  using  three  different  gas  models;  perfect 
gas,  equilibrium  air,  and  nonequilibrium  air.  The  Roe  scheme 
obtained  a  converged  steady-state  solution  in  about  3000 
iterations  using  an  81  x  41  grid,  and  a  Courant  number  of 
0.9.  On  the  Cray-YMP  this  took  only  about  2  minutes. 
Molecular  vibration  modeling  did  not  appreciably  impact 
convergence  rate.  This  robustness  is  impressive  in  spite  of 
the  complications  associated  with  the  axisymmetric  equations 
of  motion.  The  scheme's  utility  for  producing  reasonable 
preliminary  solutions  is  certainly  supported  by  the  results  of 
the  present  effort. 

The  perfect  gas  model  predicts  the  proper  qualitative 
trend  in  shock  standoff  distance,  for  inviscid  flow  over  a 
sphere.  By  incorporating  an  equilibrium  vibration  model,  both 
shock  standoff  distance  and  temperature  decrease  by  about  10- 
20  percent. 

Changing  thermal  states  near  the  shock  were  determined 
using  nonequilibrium  modeling.  By  calculating  these  rapidly 
changing  thermal  states,  the  nonequilibrium  model  predicted 
the  highest  temperature,  and  lowest  density  near  the  shock. 


A.  Solution  Sensitivity  to  Axisvmmetric  Equations  of  Motion 

Since  the  2-D  investigation  did  not  encounter  any 
spurious  oscillations,  this  anomaly  is  unique  to  the 
axisymmetric  problem.  For  the  axisymmetric  configuration 
there  is  a  geometric  singularity  and  the  Jacobian  vanishes 
along  the  axis  of  symmetry.  The  oscillatory  behavior  did  not 
affect  the  final  steady-state  solution.  In  fact,  the  noise 
near  the  symmetry  line  could  be  filtered  out  using  post¬ 
processing.  Therefore  the  behavior  was  not  considered  to  be 
a  serious  determinant  to  this  analysis. 

Nonphysical  recirculations  were  avoided  by  constructing 
a  grid  with  azimuthal  spacing  obtained  through  trial  and 
error.  No  concrete  justification  for  this  can  be  offered. 
More  investigation  into  this  phenomenon  is  warranted. 

B.  Accuracy  Versus  Computational  Complexity 

The  Millikan  and  White  model  accurately  predicts 
relaxation  rates  for  the  flight  regime  analyzed.  The  Park 
model  has  several  coefficients  which  could  be  modified  to 
bring  it  into  closer  agreement  with  the  experimental  data  used 
in  this  study.  Due  to  its  simplicity  the  Millikan  and  White 
model  is  preferred.  However,  if  the  models  are  used  over  a 
greater  variation  of  Mach  numbers,  the  Park  model  may  prove 
more  accurate  at  the  higher  Mach  numbers  (5:7,32:234). 


since  a  flow  in  thermal  nonequilibrium  develops  a  "hot 
spot"  just  downstream  of  the  shock,  for  moderate  Mach  numbers, 
equilibrium  may  not  be  reached  in  the  stagnation  region. 
Moreover,  moderate  Mach  numbers  pose  a  unique  problem  for  the 
"two-temperature"  model.  At  lower  temperatures  (i.e.  lower 
Mach  numbers) ,  the  vibrational  temperatures  for  oxygen  and 
nitrogen  at  the  body  surface  are  significantly  different.  At 
higher  temperatures  (i.e.  higher  Mach  numbers),  equilibrium  is 
attained  at  the  body  surface,  but  magnitudes  for  nitrogen  and 
oxygen  differed  by  15  percent  within  the  relaxation  zone. 
These  two  factors  may  impact  heat  transfer  calculations  and 
calculations  involving  chemical  reactions  at  the  onset  of 
dissociation. 

Equilibrium  calculations  neglect  the  relaxation  of 
temperature  and  density  within  the  shock  layer.  However,  ;.ear 
the  body  surface  equilibrium  and  nonequilibrium  calculations 
differed  by  only  about  3-5  percent  above  Mach  6.  Therefore, 
in  the  absence  of  dissociation,  equilibrium  can  be  used  to 
accurately  predict  inviscid  surface  conditions. 

Unfortunately,  the  increase  in  CPU  time  required  to 
compute  an  equilibrium  solution  is  intricately  linked  with  the 
Roe  averaging.  This  makes  it  less  desirable  than  the 
nonequilibrium  solution. 


C.  Recommendations  for  Follow-on  Study 


1.  Extension  of  the  thermal  nonequilibrium  analysis  to 
incorporate  a  second-order-accurate  TVD  scheme  is  an 
appropriate  next  step.  Here  the  impact  of  sophisticated 
dissipation  functions  can  be  assessed. 

2.  Further  analysis  of  the  grid/boundary  condition 
interaction  may  shed  further  light  on  the  development  of 
spurious  behavior  near  the  symmetry  line.  Particularly,  the 
limiting  form  of  the  equations  of  motion  should  be  more  fully 
analyzed. 

3.  More  detailed  analysis  using  a  multi-temperature 
model  with  the  Roe  scheme  would  help  determine  if  the  two- 
temperature  model  is  sufficient  at  Mach  numbers  between  5-10. 
The  multi-temperature  model  is  outlined  in  Reference  5,  and 
carried  out  for  Mach  25-28  flows  using  a  Gauss-Seidel  line 
relaxation  technique. 

4.  Extension  of  first-order,  inviscid  analysis  to  first- 
order  viscous  analysis.  Spurious  oscillations  may  go  away  by 
introducing  "natural"  dissipation  into  the  problem. 
Additionally,  the  presence  of  a  boundary-layer  should  produce 
more  accurate  shock  standoff  distances.  Finally,  computed 
heat  transfer  rates  can  be  compared  to  experiment  to  further 
assess  the  impact  of  each  molecular  vibration  model. 
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Appendix  A;  Derivations  of  Governing  Ecmations 


Perfect  Gas/Equilibrium 

The  equations  of  motion  for  2-D  or  axisymmetric  flow  can 
be  written  in  conservation  form  as: 


puy“ 

(P+pu^)  y“ 
puvy® 

{P+E^) 

pvy* 
puvy“ 

(P+v^)y‘ 
iP+E,)  vy“ 

For  2-D  flow  a  is  0,  and  for  axisymmetric  flow  a  is  1.  To 
apply  these  equations  in  a  general  coordinate  system 
ot  5  and  T]  ,  we  use  the  chain-rule: 

^x'  ^y>  ’lx'  ’ly  familiar  metrics  of  the  transformation. 

By  multiplying  through  by  the  Jacobian  of  the  inverse 
transformation,  J  ,  the  equation  can  be  recast  into  the 
"strong  conservation  form." 


where, 


U  = 


S  = 


py* 

puy“ 

pvy“ 

P.y* 

0 

0 

Pa 

0 


E  = 


F  = 


First, 


JU,  +  -h  nA  +  \yF^  +  Tl^„)  =  JS 

Assuming  that  the  grid  is  fixed  in  time,  JU^  =  (JU)  ^ 


= 

[^A( 

= 

= 

[j?A« 

-  F(J^Pj 

= 

[^TlAn 

Note  that, 

=  y^5  -  yj^  =  0 


In  a  similar  fashion. 


+  (^Vr,  =  0 

Thus, 

{JU)  ^  = 

(Jo),  +  +  [J(ti^  +  Tlj^)]^  = 


Finally,  if  we  let  U  =  JU  and. 


Also, 


JS 

JS 


E  =  J(i^  +  iyF) 


F  =  j{r\^  +  ti^y) 


Now  the  governing  equations  can  be  written  in  the  form 
suggested  by  Viviand  and  Vinokur: 


0 


I 


There  are  two  drawbacks  to  this  form.  One  drawback  of  this 
form  is  that  new  Jacobian  matrices  (and  consequently 
R  and  A  )  must  be  developed,  since  U  is  now  the  vector  of 
dependent  variables  and  not  U.  The  more  general  Jacobian 
matrices  are  given  by: 


A  =  and  B 

dU 


dU 


Since  E  and  F  are  most  conveniently  expressed  in  terms  of  the 
elements  of  U  ,  A  and  B  are  also  determined  by  using  the 
chain  rule: 

=  fuM-  wi 


Similarly 

B  =  +  11^ 

The  second  drawback  is  that  a  grid  with  varying  metrics  will 
not  conserve  the  freestream.  Consider  E,F  constant.  This 
implies  E  and  F  are  constant,  if  and  only  if  y,,,  y^,  are 

constant. 
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The  problem  can  be  avoided  by  simply  writing  the 


governing  equations  in  the  following  form  (38:2): 


The  solution  can  be  obtained  by  two  successive  sweeps,  one  in 
the  $  direction  the  other  in  the  i\  direction.  In  discretized 
form,  the  equation  for  the  t|  -sweep  becomes: 


-  Fh)]  + 


where  the  dissipation  term  is  added  for  stability.  For  the 
Roe  scheme  the  dissipation  term  is  defined  as: 


and 


f)" 

^i*X/2,j 


=  R/  |A,I  Ri‘ 


=  R<  ia;i  r/ 


The  variables  used  in  R^,  |A^|,  are  found  by  using  the  Roe- 
average  at  the  interfaces  (i+l/2,j)  and  (i-l/2,j).  The 
matrix,  |A^|,  is  a  diagonal  matrix  of  eigenvalues.  The 
derivation  of  the  linearized  matrices  is  performed  in  Appendix 
B.  Since  the  vector  of  unknowns  is  still  U  ,  the  definition 
of  A  remains  the  same  as  discussed  above  for  the  strong 


conservation  form. 


Noneaui librium 


Here  an  equation  is  included  for  each  species  which  is 
not  in  equilibrium  (9:7).  The  conseirvation  form  is  again 
written  as: 

U,  ^  Fy  =  S 


However,  now  for  a  system  with  say  two  species,  diatomic 
oxygen  and  nitrogen  for  example. 


• 

po,y“ 

Po,^^y“ 

(P+pu^)y“ 

u  = 

p  vy* 

E  = 

puvy* 

• 

(P+E^)  uy'^ 

E^^uy^ 

• 

0 

0 

Po,vy“ 

0 

puvy* 

s  = 

Pa 

F  = 

(p+v2)y« 

• 

0 

(P+PJ  vy“ 

II 

K  y“ 

where  is  the  vibrational  energy  of  species  s,  and  K,  is 
the  Landau-Teller  form  for  the  energy  exchange  between 
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translational  and  vibrational  energy  inodes,  also  for  species 
s.  The  expansion  for  a  general  coordinate  system  then 
precedes  as  above. 

The  derivation  of  the  linearized  matrices  used  in  the 
Roe-averaged  dissipation  term,  followed  the  work  of  Liu  and 
Vinokur  (25) .  Details  are  again  provided  in  Appendix  B. 


Appendix  Bs  Derivation  of  Flux  Jacobians,  Eigenvalues,  and 
Eigenvectors 

To  study  high-temperature  effects,  Roe's  approximate 
Riemann  solver  was  selected  because  it  is  both  simple  and 
robust.  The  "Roe-averaged"  state  was  derived  by  Roe  for  a 
perfect  gas  (35) ,  and  generalized  for  both  equilibrium  and 
nonequilibrium  flow  by  Vinokur  and  Liu  (24,42). 

Perfect  GasXEouilibrium 

Consistent  with  the  derivation  by  Liu  and  Vinokur  the 
flux  Jacobians  for  perfect  gas  or  equilibrium  analysis  are  as 
follows : 

1  0  O' 

(2-k)u  -kv  k 

V  u  0 

H-ku^  -xuv  H+k)  u 

0  1  o' 

V  u  0 

-KU  (2-K)  V  K 

-KUV  H-KU^  (1+K)v 
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A  = 


dE 

dU 


K^-u^ 


-uv 


B  = 


dF 

dU 


-uv 


K^-v^ 


(K^-H)  u 


r 


where 


iq  =  -|  iu^+v^)  +  X 


K 

X 


dP 

die) 


•  dP 

.dp 


Is 


e  =  pe 


For  perfect  gas  flows,  k  is  equal  to  y-l  and  is  constant. 
Also  for  perfect  gas  flows  X  is  exactly  zero. 

I  used  Mathematica  to  derive  the  eigenvectors  for  a 
general  coordinate  system  by  letting  A  =  +  ^yB  ; 

[  1  1  1  0  1 


u-k.^c 

V-K^C 


V  u+k^c  -k^ 

V  v+k2  c  k^ 


iH-k^uc-k2Vc) 


{-—+JLJJL.)  iH+k2UC+k2Vc) 

K  2  12 


ik^v-k2u) 


/  k,u  A,v^ 

(  k,  \ 

K  1  1  A  2 

—  h  77—  ^ 

^  c  c 

^  c 

1  2  J 

1  2 

l-i?i  b2U 


-b2V- 


A,  \ 


1  /-I 


k2U-k^v 


/  .  \ 
-bjU+ — - 

^  c 


-A, 


(  ,  A,  N 

-h  v+-i- 
“  c 

[h] 

1  2  J 

[  2  j 

A^  0 


77 


where 


The  eigenvalues  for  the  above  system  are 


Ax 


U^-K^C  0  0  0 

0  0  0 
0  0  U^-^K^c  0 

0  0  0  c/^ 
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» 


» 


» 


Nonecfuilibrium 

For  nonequi librium  flow  with  two  species  the  system  has 
a  rank  of  seven.  Again,  the  equation  was  transformed  to  a 
general  coordinate  system  then  the  eigenvalues  and 
eigenvectors  of  the  flux  Jacobians  were  determined  with 
Mathematica.  The  procedure,  while  involving  significant 
analytical  complexity,  is  straight  forward.  Extension  for 
mixtures  with  additional  species  is  also  straight  forward. 
These  matrices  are  provided  as  follows: 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0  0  0 

0  0  0 

0  0 

0  0 

0  0 

0  0  0 

0  0  0 


0  0 

0  0 

0  0 

0  0 

0  0 

U^-K^c  0 


0 
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(ItjV-^jU)  1  1  — 2 — ^  iH-k^uc-kjVc)  ^H+k^uc+kj^vc) 

0  0  1  0  0  e,Ci  ©iC, 

0  1  0  0  0  SjCj  ©jCj 
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iCjU-JCjV 
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l+SjC 

-6,c,b, 
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u(e^Cib,) 

vie^Cyb^) 
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eiCjj 

-Cibj 
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Wcjbj) 

-cA 

c^by 

C’r-t'. 

b.  -h  -h 
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2  ~  2  2 
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where , 
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X.U  XjV' 

i>.  +  -1-  +  — i— 
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-m-4] 
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^  c 

k,v^ 
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'  k^u  k^v' 
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2 
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2 

2 

The  present  study  was  concerned  with  only  single  species 
systems.  Therefore,  the  first  and  second  rows  and  columns 
of  |A^|  are  eliminated;  the  second  and  seventh  rows  and  second 
and  fifth  columns  in  %  are  eliminated;  the  second  and  fifth 
rows  and  second  and  seventh  columns  of  are  eliminated. 
This  reduces  the  rank  of  the  system  from  seven  to  five. 
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The  effects  of  thermal  nonequilibrium  on  flows  about 
blunt-bodies  have  not  been  studied  numerically  in  isolation  from 
chemical  reactions.  Typically,  air  is  modeled  as  a  perfect  gas  or  as 
a  chemically  reacting  mixture.  In  the  former  case,  significant  errors 
result  at  Mach  numbers  exceeding  about  5.  However,  below  Mach  numbers 
around  8,  the  effects  of  chemical  reactions  are  negligible.  This 
study  analyzes  flow  about  a  simple  axisymmetric  blunt-body  at  moderate 
hypersonic  speeds  (Mach  numbers  between  5-8) .  A  first-order-accurate 
Roe  scheme  was  used  to  compute  solutions  of  the  Euler  equations 
assuming  different  gas  models:  perfect  gas,  thermal  equilibrium,  and 
thermal  nonequilibrium.  Equilibrium  and  nonequilibrium  shock  standoff 
distances  deviated  10-15  percent  from  perfect  gas  calculations. 
Nonequilibrium  solutions  produced  striking  temperature  and  density 
gradients  in  th^  stagnation  region  with  a  unique  "hot  spot"  near  the 
shock .  '  ^  -  r  ■ 
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